#install.packages("rJava")
# Paquete para acceder datos en GBIF
#install.packages("rgbif")
# Paquete para acceder datos geoespaciales
#install.packages("geodata")
# Paquete para mapas interactivos
#install.packages("leaflet")
# Paquete para modelado de distribución de especies
#install.packages("dismo")
# install.packages("RColorBrewer")Tarea 3 Modelo de nicho de la especie Pharomachrus mocinno (Quetzal)
Instalación de paquetes
Esto se hace en la terminal
Carga de paquetes
# Colección de paquetes de Tidyverse
library(tidyverse)Warning: package 'tidyverse' was built under R version 4.3.3
Warning: package 'ggplot2' was built under R version 4.3.3
Warning: package 'tidyr' was built under R version 4.3.3
Warning: package 'readr' was built under R version 4.3.3
Warning: package 'purrr' was built under R version 4.3.3
Warning: package 'stringr' was built under R version 4.3.3
Warning: package 'forcats' was built under R version 4.3.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Estilos para ggplot2
library(ggthemes)Warning: package 'ggthemes' was built under R version 4.3.3
# Paletas de colores de RColorBrewer
library(RColorBrewer)
# Paletas de colores de viridis
library(viridisLite)Warning: package 'viridisLite' was built under R version 4.3.3
# Gráficos interactivos
library(plotly)Warning: package 'plotly' was built under R version 4.3.3
Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
# Manejo de datos vectoriales
library(sf)Warning: package 'sf' was built under R version 4.3.3
Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
# Manejo de datos raster
library(terra)Warning: package 'terra' was built under R version 4.3.3
terra 1.7.83
Attaching package: 'terra'
The following object is masked from 'package:tidyr':
extract
# Manejo de datos raster
library(raster)Warning: package 'raster' was built under R version 4.3.3
Loading required package: sp
Warning: package 'sp' was built under R version 4.3.3
Attaching package: 'raster'
The following object is masked from 'package:plotly':
select
The following object is masked from 'package:dplyr':
select
# Mapas interactivos
library(leaflet)Warning: package 'leaflet' was built under R version 4.3.3
# Acceso a datos en GBIF
library(rgbif)Warning: package 'rgbif' was built under R version 4.3.3
# Datos geoespaciales
library(geodata)Warning: package 'geodata' was built under R version 4.3.3
# Modelado de distribución de especies
library(dismo)Warning: package 'dismo' was built under R version 4.3.3
# Instalar y cargar el paquete RColorBrewer si no lo tienes
library(RColorBrewer)Selección de la especie
# Nombre de la especie
especie <- "Pharomachrus mocinno La Llave, 1832"
# Consulta a GBIF
respuesta <- occ_search(
scientificName = especie,
hasCoordinate = TRUE,
hasGeospatialIssue = FALSE,
limit = 10000
)Extracción y guardado de datos
# Extraer datos de presencia
presencia <- respuesta$data
# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')Lectura de los datos de especie
# Leer los datos de presencia de un archivo CSV
presencia <- read_csv('presencia.csv')Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
dat <- vroom(...)
problems(dat)
Rows: 10000 Columns: 116
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (77): scientificName, issues, datasetKey, publishingOrgKey, installatio...
dbl (24): key, decimalLatitude, decimalLongitude, crawlId, taxonKey, kingdo...
lgl (9): isSequenced, isInCluster, eventID, identificationRemarks, organis...
dttm (6): lastCrawled, lastParsed, dateIdentified, eventDate, modified, las...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
presencia <- st_as_sf(
presencia,
coords = c("decimalLongitude", "decimalLatitude"),
remove = FALSE, # conservar las columnas de las coordenadas
crs = 4326
)Gráfico de la distribución por países de la especie Pharomachrus mocinno
# Gráfico ggplot2
grafico_ggplot2 <-
presencia |>
st_drop_geometry() |>
ggplot(aes(x = fct_infreq(countryCode))) +
geom_bar(
aes(
text = paste0(
"Cantidad de registros de presencia: ", after_stat(count)
)
)
) +
ggtitle("Cantidad de registros de Pharomachrus mocinno por país") +
xlab("País") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: GBIF") +
theme_economist()Warning in geom_bar(aes(text = paste0("Cantidad de registros de presencia: ", :
Ignoring unknown aesthetics: text
# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')Gráfico de Distribución porcentual tipo pastel de la especie Pharomachrus mocinno por Provincia en Costa Rica
# Filtrar datos para Costa Rica
presencia_cr <- presencia |>
filter(countryCode == "CR")
# Calcular porcentajes de distribución por provincia
distribucion_provincia <- presencia_cr |>
st_drop_geometry() |>
group_by(stateProvince) |>
summarise(count = n()) |>
mutate(percentage = (count / sum(count)) * 100)
# Filtrar el dataset para eliminar las provincias San José, Guanacaste y Puntarenas
distribucion_provincia_filtrada <- distribucion_provincia |>
filter(!stateProvince %in% c("Provincia San José", "Provincia Guanacaste", "Provincia Puntarenas"),
!is.na(stateProvince))# Crear gráfico tipo pie con ggplot2
grafico_pie <-
distribucion_provincia_filtrada |>
ggplot(aes(
x = "",
y = percentage,
fill = stateProvince,
text = paste0(
"Provincia: ", stateProvince,
"<br>Porcentaje: ", round(percentage, 1), "%"
)
)) +
geom_bar(stat = "identity", width = 1, color = "white") +
coord_polar(theta = "y") +
ggtitle("Distribución de registros de Pharomachrus mocinno \n por provincia en Costa Rica") +
geom_text(
aes(label = paste0(round(percentage), "%")),
color = "black",
size = 3, # Ajustar el tamaño del texto
position = position_stack(vjust = 0.5) # Para ajustar la posición del texto en cada porción
) +
scale_fill_brewer(palette = "Paired") + # Cambiar la paleta de colores (puedes probar con otras paletas)
theme_void() +
theme(
plot.title = element_text(hjust = 0.5),
legend.title = element_blank()
)
# Mostrar gráfico
grafico_pie
Mapa de la distribución de la especie Pharomachrus mocinno
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Amanita muscaria"
) |>
addLayersControl(
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c("Registros de Bradypus variegatus"))Mapa de la distribución de la especie Pharomachrus mocinno junto con variables climáticas de precipitación y temperatura
# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())# Nombres de las variables climáticas
names(clima) [1] "wc2.1_10m_bio_1" "wc2.1_10m_bio_2" "wc2.1_10m_bio_3" "wc2.1_10m_bio_4"
[5] "wc2.1_10m_bio_5" "wc2.1_10m_bio_6" "wc2.1_10m_bio_7" "wc2.1_10m_bio_8"
[9] "wc2.1_10m_bio_9" "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
[13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
[17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"
# Definir la extensión del área de estudio
area_estudio <- ext(
min(presencia$decimalLongitude) - 5,
max(presencia$decimalLongitude) + 5,
min(presencia$decimalLatitude) - 5,
max(presencia$decimalLatitude) + 5
)# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
# palette = "inferno",
# palette = "magma",
palette = rev(brewer.pal(11, "RdYlBu")),
values(clima$wc2.1_10m_bio_1),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
# palette = "viridis",
# palette = "YlGnBu",
palette = "Blues",
values(clima$wc2.1_10m_bio_12),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
clima$wc2.1_10m_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
clima$wc2.1_10m_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación",
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLegend(
title = "Temperatura",
values = values(clima$wc2.1_10m_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(clima$wc2.1_10m_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c("Temperatura", "Precipitación", "Registros de Pharomachrus mocinno")
) |>
hideGroup("Precipitación")Creación del modelo de nichos para la especie Pharomachrus mocinno
# Crear dataframe con columnas de longitud y latitud
coordenadas_presencia <- data.frame(
decimalLongitude = presencia$decimalLongitude,
decimalLatitude = presencia$decimalLatitude
)
# Eliminar coordenadas duplicadas
coordenadas_presencia <- unique(coordenadas_presencia)# Establecer una "semilla" para garantizar que la selección aleatoria sea reproducible
set.seed(123)
# Cantidad de registros de presencia
n_presencia <- nrow(coordenadas_presencia)
# Con sample(), se selecciona aleatoriamente una proporción (ej. 0.7)
# de los índices de los datos de presencia para el conjunto de entrenamiento
indices_entrenamiento <- sample(
1:n_presencia,
size = round(0.7 * n_presencia)
)
# Crear el subconjunto de entrenamiento utilizando los índices seleccionados
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
# Crear el subconjunto de evaluación con los datos restantes
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima <- raster::stack(clima)
# Ejecutar el modelo
modelo_maxent <- maxent(x = clima, p = entrenamiento)Loading required namespace: rJava
# Aplicar el modelo entrenado a las variables climáticas
# para generar un mapa de idoneidad del hábitat
prediccion <- predict(modelo_maxent, clima)# terra::extract() extrae los valores del raster de predicción
# en las coordenadas de evaluación
# eval_pres almacena los valores de idoneidad predichos
# en los puntos de evaluación de presencia
eval_pres <- terra::extract(
prediccion,
evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)
# Generar puntos aleatorios dentro del área de estudio definida.
# Estos puntos se asumen como ausencias de la especie.
ausencias <- randomPoints(mask = clima, n = 1000)
# eval_aus almacena los valores de idoneidad predichos
# en los puntos de ausencia
eval_aus <- terra::extract(
prediccion,
ausencias
)
# Generar estadísticas de evaluación del modelo
resultado_evaluacion <- evaluate(p = eval_pres, a = eval_aus)Curva ROC y resultado del dato AUC del modelo par ala especie Pharomachrus mocinno
# Datos para graficar la curva ROC
datos_roc <- data.frame(
FPR = resultado_evaluacion@FPR,
TPR = resultado_evaluacion@TPR,
Umbral = resultado_evaluacion@t
)
# Valor AUC
auc <- resultado_evaluacion@auc
# Gráfico ggplot2
grafico_ggplot2 <-
ggplot(
datos_roc,
aes(
x = FPR,
y = TPR,
u = Umbral
)
) +
geom_line(
color = "blue",
size = 1
) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
labs(title = paste("Curva ROC (AUC =", round(auc, 3), ")"),
x = "Tasa de falsos positivos (FPR)",
y = "Tasa de verdaderos positivos (TPR)") +
theme_minimal()Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
# Gráfico plotly
ggplotly(grafico_ggplot2) |>
config(locale = 'es')Mapa interactivo de idoneidad para la probabilidad de encontrar la especie Pharomachrus mocinno
# Paleta de colores del modelo
colores_modelo <- colorNumeric(
palette = c("white", "black"),
values(prediccion),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
clima$wc2.1_10m_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
clima$wc2.1_10m_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación",
) |>
addRasterImage( # capa raster del modelo de distribución
prediccion,
colors = colores_modelo,
opacity = 0.6,
group = "Modelo de distribución",
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLegend(
title = "Temperatura",
values = values(clima$wc2.1_10m_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(clima$wc2.1_10m_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Modelo de distribución",
values = values(prediccion),
pal = colores_modelo,
position = "bottomright",
group = "Modelo de distribución"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Temperatura",
"Precipitación",
"Modelo de distribución",
"Registros de Pharomachrus mocinno"
)
) |>
hideGroup("Temperatura") |>
hideGroup("Precipitación")Mapa de idoneidad binario para la especie Pharomachrus mocinno
# Definir el umbral
umbral <- 0.5
# Crear el raster binario
prediccion_binaria <- (prediccion >= umbral) * 1
# Crear la paleta de colores para el raster binario
colores_prediccion_binaria <- colorFactor(
palette = c("transparent", "blue"), # "transparent" para las áreas no adecuadas
domain = c(0, 1),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage(
prediccion_binaria,
colors = colores_prediccion_binaria,
opacity = 0.6,
group = "Modelo de distribución binario",
) |>
addCircleMarkers(
data = presencia,
stroke = FALSE,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLegend(
title = "Modelo de distribución binario",
labels = c("Ausencia", "Presencia"),
colors = c("transparent", "blue"),
position = "bottomright",
group = "Modelo de distribución binario"
) |>
addLayersControl(
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Modelo de distribución binario",
"Registros de Pharomachrus mocinno"
)
)Warning in colors(.): Some values were outside the color scale and will be
treated as NA